function [out_xx,out_yy]=mcgh_EMforGaussian(data,binwidth,row);
global p_1 p1 p0 m1 m0 m_1 s_1 s1 s0
%use EM to estimate the parameter
    %initial value
    x=data;
    rep_x=repmat(x,1,3);

    len_x=length(x);
    alpha=[p_1 p0 p1];
    rep_alpha=repmat(alpha,len_x,1);

    sigma=[s_1 s0 s1];
    rep_sigma=repmat(sigma,len_x,1);

    mue=[m_1 m0 m1];
    rep_mue=repmat(mue,len_x,1);
    %Q_old=0;
    idx=1;
    
    %pn=n./sum(n);
    %hh=bar(xbin,n,1);

    %xd = get(hh,'Xdata'); % Gets the x-data of the bins.
    %np = get(gca,'NextPlot');    
    %set(gca,'NextPlot','add')                                 
    %xd = get(hh,'Xdata'); % Gets the x-data of the bins.
    %rangex = max(xd(:)) - min(xd(:)); % Finds the range of this data.
    %binwidth = rangex/nbins;    % Finds the width of each bin.

    xx=[];
    yy=[];		
    while 1
   
        %three gauss %alpha:pi, sigma: standard deviation, mue: mean of center
        %f_x=rep_alpha.*((2.*pi.*rep_sigma).^-0.5).*exp(-0.5.*(rep_x-rep_mue).^2./(rep_sigma.^2));
        %f_x(find(isnan(f_x)==1))=0

        f_x(:,1)=alpha(1).*mcgh_normpdf2(x,mue(1),sigma(1));
        f_x(:,2)=alpha(2).*mcgh_normpdf2(x,mue(2),sigma(2));
        f_x(:,3)=alpha(3).*mcgh_normpdf2(x,mue(3),sigma(3));
   
        f_x = row*(f_x*binwidth);   % Normalization necessary to overplot the histogram.
        f_x=f_x./row;
        sum_f_x=sum(f_x,2); 
    
        %E step
        log_f_x=log(f_x);
	%pause
        log_f_x(isinf(log_f_x)==1)=0;
    
        log_sum_f_x=log(sum_f_x);
        log_sum_f_x(isinf(log_sum_f_x)==1)=0;
        log_sum_f_x(isnan(log_sum_f_x)==1)=0;
	log_y_ij_g=log(f_x)-repmat(log_sum_f_x,1,3);
        y_ij_g=exp(log_y_ij_g);
    
        y_ij_g(find(isnan(y_ij_g)==1))=0;
        Q(idx)=sum(sum(y_ij_g.*log(rep_alpha)))+sum(sum(y_ij_g.*log_f_x));
       
	 if (idx >1 & Q(idx)<Q(idx-1)) | isnan(Q(idx))
            break;
        else
            p_1=alpha(1);
            p0=alpha(2);
            p1=alpha(3);
            s_1=sigma(1);
            s0=sigma(2);
	    s1=sigma(3);
            m_1=mue(1);
            m0=mue(2);
            m1=mue(3);
	    out_xx=xx;
	    out_yy=yy;	
        end
    
        %plot fit the gaussina
        %clf
        %hh=bar(xbin,n,1);

        mstd=std(x);
        %added june2004
        [uniq_x uidx]=unique(x);
        xx=min(min(uniq_x)):0.1*mstd:max(max(uniq_x));
        yy=interp1(uniq_x,sum_f_x(uidx),xx,'nearest');
        %end added Jun2 2004
        idx=idx+1;
    
    
        %M step
        alpha_j_g1=sum(y_ij_g,1)/len_x;
        alpha=alpha_j_g1;
    
        mue_j_g1=1./len_x./alpha_j_g1.*sum(rep_x.*y_ij_g,1);
        mue=mue_j_g1;
    
        sigma_j_g1=(1/len_x./alpha_j_g1.*sum(y_ij_g.*(rep_x-repmat(mue_j_g1,len_x,1)).^2,1)).^0.5;
        sigma=sigma_j_g1;
    
        rep_alpha=repmat(alpha_j_g1,len_x,1);
        rep_sigma=repmat(sigma_j_g1,len_x,1);
        rep_mue=repmat(mue_j_g1,len_x,1);
        % pause
    end %end while
    %p_1, p1, p0,
    %m1, m0, m_1, 
    %s_1, s1, s0
